library(survey)
library(ggplot2)

## run "PCM_cleaning.R"

## Create weighted survey design object
milconf_design <-
  svydesign(
    id = ~ 1,
    weights = ~ weight,
    data = mil_conf.df
  )

#### CONFIDENCE: Q9 BY CONDITION ####

Q9_diff <- matrix(NA, ncol = 3, nrow = 7)

## Create difference in means for each condition
diff_Q9_2 <- svyglm(Q9_m ~ ASSIGNMENT_A_2, milconf_design)
Q9_diff[1,1] <- diff_Q9_2$coefficients[2]
Q9_diff[1,2] <- confint(diff_Q9_2)[2,1]
Q9_diff[1,3] <- confint(diff_Q9_2)[2,2]

diff_Q9_3 <- svyglm(Q9_m ~ ASSIGNMENT_A_3, milconf_design)
Q9_diff[2,1] <- diff_Q9_3$coefficients[2]
Q9_diff[2,2] <- confint(diff_Q9_3)[2,1]
Q9_diff[2,3] <- confint(diff_Q9_3)[2,2]

diff_Q9_4 <- svyglm(Q9_m ~ ASSIGNMENT_A_4, milconf_design)
Q9_diff[3,1] <- diff_Q9_4$coefficients[2]
Q9_diff[3,2] <- confint(diff_Q9_4)[2,1]
Q9_diff[3,3] <- confint(diff_Q9_4)[2,2]

diff_Q9_5 <- svyglm(Q9_m ~ ASSIGNMENT_A_5, milconf_design)
Q9_diff[4,1] <- diff_Q9_5$coefficients[2]
Q9_diff[4,2] <- confint(diff_Q9_5)[2,1]
Q9_diff[4,3] <- confint(diff_Q9_5)[2,2]

diff_Q9_6 <- svyglm(Q9_m ~ ASSIGNMENT_A_6, milconf_design)
Q9_diff[5,1] <- diff_Q9_6$coefficients[2]
Q9_diff[5,2] <- confint(diff_Q9_6)[2,1]
Q9_diff[5,3] <- confint(diff_Q9_6)[2,2]

diff_Q9_7 <- svyglm(Q9_m ~ ASSIGNMENT_A_7, milconf_design)
Q9_diff[6,1] <- diff_Q9_7$coefficients[2]
Q9_diff[6,2] <- confint(diff_Q9_7)[2,1]
Q9_diff[6,3] <- confint(diff_Q9_7)[2,2]

diff_Q9_8 <- svyglm(Q9_m ~ ASSIGNMENT_A_8, milconf_design)
Q9_diff[7,1] <- diff_Q9_8$coefficients[2]
Q9_diff[7,2] <- confint(diff_Q9_8)[2,1]
Q9_diff[7,3] <- confint(diff_Q9_8)[2,2]

## Plot
Q9_plot.df <- data.frame(change = Q9_diff[,1],
                         ci_low = Q9_diff[,2],
                         ci_hi = Q9_diff[,3],
                         cond = c(1:7))

Q9_plot <- ggplot(Q9_plot.df) + 
  geom_bar(aes(x = cond, y = change), stat = "identity", width = 0.75, fill = "gray60") +
  geom_errorbar(aes(x = cond, ymin = ci_low, ymax = ci_hi), width = 0.4)
Q9_plot <- Q9_plot + xlab("Treatment Condition") + 
  scale_x_continuous(breaks = c(1:7),
                     labels = c("2", "3", "4", "5", "6", "7", "8")) +
  ylab("Difference in Means from Control Condition") + ggtitle("Q9 by Treatment Condition")
Q9_plot <- Q9_plot + geom_hline(yintercept = 0, colour = "black", lty = 1)
Q9_plot + theme_bw()

#### CONFIDENCE: Q9 BY CONDITION AND PARTY ####

Q9_diff_p <- matrix(NA, ncol = 9, nrow = 7)

## Create difference in means, CIs for republicans
diff_Q9_2_r <- svyglm(Q9_m_rep ~ ASSIGNMENT_A_2, milconf_design)
Q9_diff_p[1,1] <- diff_Q9_2_r$coefficients[2]
Q9_diff_p[1,2] <- confint(diff_Q9_2_r)[2,1]
Q9_diff_p[1,3] <- confint(diff_Q9_2_r)[2,2]

diff_Q9_3_r <- svyglm(Q9_m_rep ~ ASSIGNMENT_A_3, milconf_design)
Q9_diff_p[2,1] <- diff_Q9_3_r$coefficients[2]
Q9_diff_p[2,2] <- confint(diff_Q9_3_r)[2,1]
Q9_diff_p[2,3] <- confint(diff_Q9_3_r)[2,2]

diff_Q9_4_r <- svyglm(Q9_m_rep ~ ASSIGNMENT_A_4, milconf_design)
Q9_diff_p[3,1] <- diff_Q9_4_r$coefficients[2]
Q9_diff_p[3,2] <- confint(diff_Q9_4_r)[2,1]
Q9_diff_p[3,3] <- confint(diff_Q9_4_r)[2,2]

diff_Q9_5_r <- svyglm(Q9_m_rep ~ ASSIGNMENT_A_5, milconf_design)
Q9_diff_p[4,1] <- diff_Q9_5_r$coefficients[2]
Q9_diff_p[4,2] <- confint(diff_Q9_5_r)[2,1]
Q9_diff_p[4,3] <- confint(diff_Q9_5_r)[2,2]

diff_Q9_6_r <- svyglm(Q9_m_rep ~ ASSIGNMENT_A_6, milconf_design)
Q9_diff_p[5,1] <- diff_Q9_6_r$coefficients[2]
Q9_diff_p[5,2] <- confint(diff_Q9_6_r)[2,1]
Q9_diff_p[5,3] <- confint(diff_Q9_6_r)[2,2]

diff_Q9_7_r <- svyglm(Q9_m_rep ~ ASSIGNMENT_A_7, milconf_design)
Q9_diff_p[6,1] <- diff_Q9_7_r$coefficients[2]
Q9_diff_p[6,2] <- confint(diff_Q9_7_r)[2,1]
Q9_diff_p[6,3] <- confint(diff_Q9_7_r)[2,2]

diff_Q9_8_r <- svyglm(Q9_m_rep ~ ASSIGNMENT_A_8, milconf_design)
Q9_diff_p[7,1] <- diff_Q9_8_r$coefficients[2]
Q9_diff_p[7,2] <- confint(diff_Q9_8_r)[2,1]
Q9_diff_p[7,3] <- confint(diff_Q9_8_r)[2,2]

## Create difference in means, CIs for independents
diff_Q9_2_i <- svyglm(Q9_m_ind ~ ASSIGNMENT_A_2, milconf_design)
Q9_diff_p[1,4] <- diff_Q9_2_i$coefficients[2]
Q9_diff_p[1,5] <- confint(diff_Q9_2_i)[2,1]
Q9_diff_p[1,6] <- confint(diff_Q9_2_i)[2,2]

diff_Q9_3_i <- svyglm(Q9_m_ind ~ ASSIGNMENT_A_3, milconf_design)
Q9_diff_p[2,4] <- diff_Q9_3_i$coefficients[2]
Q9_diff_p[2,5] <- confint(diff_Q9_3_i)[2,1]
Q9_diff_p[2,6] <- confint(diff_Q9_3_i)[2,2]

diff_Q9_4_i <- svyglm(Q9_m_ind ~ ASSIGNMENT_A_4, milconf_design)
Q9_diff_p[3,4] <- diff_Q9_4_i$coefficients[2]
Q9_diff_p[3,5] <- confint(diff_Q9_4_i)[2,1]
Q9_diff_p[3,6] <- confint(diff_Q9_4_i)[2,2]

diff_Q9_5_i <- svyglm(Q9_m_ind ~ ASSIGNMENT_A_5, milconf_design)
Q9_diff_p[4,4] <- diff_Q9_5_i$coefficients[2]
Q9_diff_p[4,5] <- confint(diff_Q9_5_i)[2,1]
Q9_diff_p[4,6] <- confint(diff_Q9_5_i)[2,2]

diff_Q9_6_i <- svyglm(Q9_m_ind ~ ASSIGNMENT_A_6, milconf_design)
Q9_diff_p[5,4] <- diff_Q9_6_i$coefficients[2]
Q9_diff_p[5,5] <- confint(diff_Q9_6_i)[2,1]
Q9_diff_p[5,6] <- confint(diff_Q9_6_i)[2,2]

diff_Q9_7_i <- svyglm(Q9_m_ind ~ ASSIGNMENT_A_7, milconf_design)
Q9_diff_p[6,4] <- diff_Q9_7_i$coefficients[2]
Q9_diff_p[6,5] <- confint(diff_Q9_7_i)[2,1]
Q9_diff_p[6,6] <- confint(diff_Q9_7_i)[2,2]

diff_Q9_8_i <- svyglm(Q9_m_ind ~ ASSIGNMENT_A_8, milconf_design)
Q9_diff_p[7,4] <- diff_Q9_8_i$coefficients[2]
Q9_diff_p[7,5] <- confint(diff_Q9_8_i)[2,1]
Q9_diff_p[7,6] <- confint(diff_Q9_8_i)[2,2]

## Create difference in means, CIs for democrats
diff_Q9_2_d <- svyglm(Q9_m_dem ~ ASSIGNMENT_A_2, milconf_design)
Q9_diff_p[1,7] <- diff_Q9_2_d$coefficients[2]
Q9_diff_p[1,8] <- confint(diff_Q9_2_d)[2,1]
Q9_diff_p[1,9] <- confint(diff_Q9_2_d)[2,2]

diff_Q9_3_d <- svyglm(Q9_m_dem ~ ASSIGNMENT_A_3, milconf_design)
Q9_diff_p[2,7] <- diff_Q9_3_d$coefficients[2]
Q9_diff_p[2,8] <- confint(diff_Q9_3_d)[2,1]
Q9_diff_p[2,9] <- confint(diff_Q9_3_d)[2,2]

diff_Q9_4_d <- svyglm(Q9_m_dem ~ ASSIGNMENT_A_4, milconf_design)
Q9_diff_p[3,7] <- diff_Q9_4_d$coefficients[2]
Q9_diff_p[3,8] <- confint(diff_Q9_4_d)[2,1]
Q9_diff_p[3,9] <- confint(diff_Q9_4_d)[2,2]

diff_Q9_5_d <- svyglm(Q9_m_dem ~ ASSIGNMENT_A_5, milconf_design)
Q9_diff_p[4,7] <- diff_Q9_5_d$coefficients[2]
Q9_diff_p[4,8] <- confint(diff_Q9_5_d)[2,1]
Q9_diff_p[4,9] <- confint(diff_Q9_5_d)[2,2]

diff_Q9_6_d <- svyglm(Q9_m_dem ~ ASSIGNMENT_A_6, milconf_design)
Q9_diff_p[5,7] <- diff_Q9_6_d$coefficients[2]
Q9_diff_p[5,8] <- confint(diff_Q9_6_d)[2,1]
Q9_diff_p[5,9] <- confint(diff_Q9_6_d)[2,2]

diff_Q9_7_d <- svyglm(Q9_m_dem ~ ASSIGNMENT_A_7, milconf_design)
Q9_diff_p[6,7] <- diff_Q9_7_d$coefficients[2]
Q9_diff_p[6,8] <- confint(diff_Q9_7_d)[2,1]
Q9_diff_p[6,9] <- confint(diff_Q9_7_d)[2,2]

diff_Q9_8_d <- svyglm(Q9_m_dem ~ ASSIGNMENT_A_8, milconf_design)
Q9_diff_p[7,7] <- diff_Q9_8_d$coefficients[2]
Q9_diff_p[7,8] <- confint(diff_Q9_8_d)[2,1]
Q9_diff_p[7,9] <- confint(diff_Q9_8_d)[2,2]

## Plot
Q9_plot_p.df <- data.frame(rep = Q9_diff_p[,1],
                           rep_lo = Q9_diff_p[,2],
                           rep_hi = Q9_diff_p[,3],
                           ind = Q9_diff_p[,4],
                           ind_lo = Q9_diff_p[,5],
                           ind_hi = Q9_diff_p[,6],
                           dem = Q9_diff_p[,7],
                           dem_lo = Q9_diff_p[,8],
                           dem_hi = Q9_diff_p[,9],
                           rcond = c(1,5,9,13,17,21,25),
                           icond = c(2,6,10,14,18,22,26),
                           dcond = c(3,7,11,15,19,23,27))

Q9_plot_p <- ggplot(Q9_plot_p.df) + 
  geom_bar(aes(x = rcond, y = rep), stat = "identity", width = 0.8, fill = "indianred1") +
  geom_errorbar(aes(x = rcond, ymin = rep_lo, ymax = rep_hi), width = 0.2) +
  geom_bar(aes(x = icond, y = ind), stat = "identity", width = 0.8, fill = "gray 83") +
  geom_errorbar(aes(x = icond, ymin = ind_lo, ymax = ind_hi), width = 0.2) +
  geom_bar(aes(x = dcond, y = dem), stat = "identity", width = 0.8, fill = "steelblue2") +
  geom_errorbar(aes(x = dcond, ymin = dem_lo, ymax = dem_hi), width = 0.2)
Q9_plot_p <- Q9_plot_p + xlab("Treatment Condition") + 
  scale_x_continuous(breaks = c(2, 6, 10, 14, 18, 22, 26),
                     labels = c("2", "3", "4", "5", "6", "7", "8")) +
  ylab("Difference in Means from Control Condition")
Q9_plot_p <- Q9_plot_p + geom_hline(yintercept = 0, colour = "black", lty = 1) +
  ggtitle("Q9 by Condition and Party ID")
Q9_plot_p + theme_bw()